Nonequilibrium dynamics of a stochastic model of 
anomalous heat transport: numerical analysis 
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+2 ! Abstract. Wc study heat transport in a chain of harmonic oscillators with random 

^ ■ clastic collisions between nearest-neighbours. The equations of motion of the 

. covariance matrix are numerically solved for free and fixed boundary conditions. In 

the thermodynamic limit, the shape of the temperature profile and the value of the 
stationary heat flux depend on the choice of boundary conditions. For free boundary 
conditions, they also depend on the coupling strength with the heat baths. Moreover, 
Q ' we find a strong violation of local equilibrium at the chain edges that determine two 

O . boundary layers of size y/~N (where N is the chain length), that are characterized 

by a different scaling behaviour from the bulk. Finally, we investigate the relaxation 
towards the stationary state, finding two long time scales: the first corresponds to the 
relaxation of the hydrodynamic modes; the second is a manifestation of the finitcness 



> 

(■ — ' of the system 



1. Introduction 

The problem of heat transport in chains of oscillators is one of the most relevant 
testing grounds to understand the behaviour of statistical systems steadily kept out 
of equilibrium. In the last decade, numerical simulations and analytic arguments have 
contributed to clarify the behaviour of such systems in the thermodynamic limit (see 
review papers [TJ [2} [3] and references therein). However, there is still a number of 
open questions such as the role of Boundary Conditions (BC in the following) and 
the convergence towards the stationary state. In spite of the continuous increase of 
computer performances, direct numerical simulations are still not so effective as to 
provide reliable data on sufficiently large systems. In this respect, stochastic models 
like the one introduced in [4] prove very helpful. In this paper we consider a version 
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of such models already analyzed in [5J [§]• The model consists in a chain of N coupled 
harmonic oscillators in interaction (at the boundaries) with two stochastic heat baths 
at different temperatures. In addition, the oscillators are subject to stochastic collisions 
that exchange the momenta of randomly chosen pairs of neighbouring oscillators, so that 
both energy and momentum are conserved. In a sense, collisions simulate the presence 
of nonlinear terms, as they contribute to ensuring ergodicity of an otherwise integrable 
model. In fact, it has been observed that this model closely reproduces the behaviour 
of standard nonlinear systems such as an FPU-/9 chain, starting from an anomalous 
(diverging) heat conductivity [5J . On the other hand, being the collision rule a perfectly 
linear process, the evolution equations for ensemble averages of the relevant observables 
can be written in an exact form and thereby solved numerically, without having to deal 
with the statistical fluctuations that affect finite samples. As a result, we have found 
that the invariant measure can be effectively approximated by the product of Gaussian 
distributions aligned along the eigendirections of the covariance matrix [5J. Moreover, 
in [6], we investigated the continuum limit (which corresponds to the large- iV limit) of 
the covariance matrix, deriving suitable partial differential equations for the stationary 
state, in the case of fixed BC. As a result, we have obtained explicit formulae for the 
temperature profile and the energy current. Remarkably, this is the first example of an 
analytic expression for the temperature profile in a system characterized by anomalous 
heat transport. 

In [7] we go beyond, by extending the continuum limit to include the time 
dependence of the covariance matrix. The reader is thus referred to [7] for a more 
detailed introduction and the corresponding bibliography. The aim of this paper is 
to complement the analysis contained in [Tj with accurate numerical studies of finite 
samples with the goal of clarifying those issues that are too difficult to be worked out 
analytically. We start by numerically computing the stationary covariance for both 
free and fixed BC. This helps to shed some light on the nontrivial role played by 
BC, whenever heat transport exhibits an anomalous behaviour. In the presence of 
normal transport, one expects that BC affect only a finite boundary layer so that, in 
the thermodynamic limit the leading term of the heat flux is independent of BC. On the 
other hand, it is known that in disordered chains of linear oscillators, the same system 
may even behave as a thermal superconductor or as an insulator, by simply switching 
from free to fixed BC [8]. In generic nonlinear chains, numerical simulations suggest 
that the heat flux scales in the same way, independently of BC. However, the careful 
simulations performed in [9] revealed that in the FPU-/5 model, the ratio between the 
heat fluxes measured for free and fixed BC does not converge to 1 for N — > oo. Here we 
show that the same behaviour occurs in our stochastic model. Actually the dependence 
on BC is even more subtle than one could have imagined: while in the case of fixed 
BC, the heat flux and the temperature profile are asymptotically independent of the 
coupling strength with the thermal baths, the same is not true for free BC. 

A second objective of this paper is the analysis of the convergence towards the 
steady state. This question, which has been hardly discussed in the literature, can be 
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straightforwardly addressed for our model, as it amounts to computing the eigenvalues of 
the evolution operator for the covariances. Moreover, we also compare the convergence 
of the average heat flux for different system sizes to show how careful direct simulations 
must be, if they have to be trusted. We find that finite-size effects associated with the 
relaxation rates of slow, i.e. long-wavelength, modes significantly modify the asymptotic 
scaling of the relaxation process. In practice, we find numerical evidence that the 
theoretical hydrodynamic scaling holds only over a finite range of time scales, although 
its duration diverges with N. 

The paper is organized as follows. In Section [2], we briefly recall the definition 
of the covariance matrix, and the coupled equations governing its evolution towards 
the stationary value. Some properties of the steady state are discussed in Section |3j 
The problem of the approach to the steady state is addressed in Section HI Finally, in 
Section [5] we summarize our main results. 

2. Equations for the covariance matrix 

In this section, we introduce the minimal notations and definitions needed to follow the 
main discussion presented in the following sections. The reader interested in a more 
detailed presentation is referred to [7J. We consider a chain of N unit-mass particles 
interacting via nearest-neighbour harmonic coupling of frequency to. The equations of 
motion are given by 

q n =Pn (1) 
p n =u 2 (5 ntN q n+ i - 2q n + <5 ni ig n _i) + 5„,i(£ + - Xqx) + S n>N (^~ - Xq N ) , 

where p n , q n are the momentum and displacement from equilibrium position of the n-th 
particle, 8ij is the Kronecker delta and 5i,j = 1 — 5i,j, and are independent Wiener 
processes with zero mean and variance 2\ksT±, where fee is the Boltzmann constant 
and A is the coupling constant. In the following free and fixed boundary conditions will 
be considered. These can be expressed in terms of the position variable q as: go = Qx, 
(In = 1n+i for free BC and q = qN+i = for fixed BC. 
We consider the covariance matrix written as 

e =(s:)- (2 » 

where the matrices y, z and v, of respective dimension (N — 1) x (N — 1), (N — 1) x iV 
and NxN are defined as 

= (AqAqj) , = (AqiPj) , V;j = (piPj) , (3) 

where (■) denotes the average over phase space probability distribution function P and 
Aqi — q{ — qi-i stand for the particle relative displacements. The variables (Aqi,pi) are 
the convenient choice to deal with: one the one hand, absolute positions are not well 
defined for free BC and on the other hand, the potential energy is expressed in terms 
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of relative differences. The only subtlety is that the domain of definition of Ag^ differs 
from that of p^. Thereby, the bulk of the system is defined as {Ag^ | i £ [2,iV]} and 
{pi | % £ [2, N — 1]}. The evolution equations for c in the bulk are 



'■j 



J ':l 



'•J 



u 2 (z i+ i 



(4) 



v,_ij + to (yy+i - yij) + 7 ( z i,j+i + z m-i 
z jti + z i+ ij - Zij) + yW it j . 



2zjj) , 



These equations follow from the deterministic equations of motion (TjQ) plus the 
contribution of the stochastic noise (7 denotes the collision rate), that is described 
by the collision matrix W, 



v j-lj-l + v i+l,j + l 

v, ;4 -i ,• + v,: - 2v. 



r i-l,j 



i.J 



V i,i+1 
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3 

J = -l 
J = l 
-J'l > 1 



(5) 



On the boundaries, several changes appear in the velocity fields. The interested reader 
can find a full description in Section 2.2 of [7j. Here we limit ourselves to show the 
contribution arising from the coupling with the heat bath, namely 



i.j 



^■lZji 



0~j,N z i,N , 



= <5j,i v *,i + S jtN v ijN + S^vxj + S i>N v Njj - 2 (T + 8 itl S jtl + T 5^ N 5 jt N ) 



3. Stationary covariance 

In this section we investigate some properties of the nonequilibrium steady state, for both 
fixed and free BC. The stationary state is obtained by considering the time-independent 
solution of equations d3J. It can be efficiently determined by exploiting the sparsity of the 
corresponding linear problem, as well as the symmetries of the unknowns (this approach 
has been followed in [5] for fixed BC). Alternatively, one can just let evolve equations (141) 
starting from any meaningful initial conditions, as the dynamics will necessarily converge 
towards the only stable stationary state (here we have adopted this latter approach also 
because we wish to study the convergence - see in the following). All numerical results 
presented in this paper have been obtained for u = 1, T + = 1.5, T_ = 0.5. This is 
by no means a limitation, as all these parameters can be easily scaled out due to the 
linear structure of the model. Accordingly, they will not be mentioned again, unless 
specifically needed for a comparison with theoretical predictions. 



3.1. The heat flux 

The first observable we have looked at is the energy flux at position % which, in terms 
of the matrices v, z is written as [8] 

Ji = -^ 2 z;+i,i+i + ^(v^ - Vj+i^+i). (6) 
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Figure 1. Free BC (A = 1): the scaled stationary flux jjN versus 1/N for 7=0.2 
(circles), 0.5 (squares), 1 (diamonds), 2 (triangles), 5 (crosses). The lines are obtained 
by fitting the data with equation ([7]). 



We have adopted the convention that a positive flux corresponds to a propagation 
towards increasing values of the spatial index i. The first term stems from the 
deterministic forces and provides for the leading (anomalous) contribution, while the 
second one accounts for energy exchanges due to collisions of nearby particles. In the 
stationary state, Jj is independent of i, i.e. Jj = J . 

In figure [TJ we show Jy/N as a function of the inverse of the system size N. The 
results refer to free BC (as we do not have analytic estimates to compare with), A = 1 and 
different values of the collision rate 7 (see the various symbols as described in the figure 
caption). In all cases there is a convincing evidence that J ~ iV -1 / 2 , similar to what 
predicted analytically in [6] for the case of fixed BC. As a consequence, the effective 
conductivity, k = JN/(T + — TJ) diverges as s/N. However, from figure [T] it is also 
evident the presence of sub-leading singular corrections which hinder the extrapolation 
of the asymptotic value. In analogy to [6j, we introduce the Ansatz, 



By using this formula to fit the data, we obtain the curves reported in figure [TJ which 
reproduce quite well the raw data. Notice that the convergence is from below for smaller 
7 values, while from above for larger collision rates. All the estimated (3 values range in 
the interval [0.88,0.95], suggesting that this parameter may be "universal". 

The extrapolated J values are plotted in figure where we can see that J ~ 7 -1 ^ 2 , 
as found for fixed BC [6] . For 7 = 5, the extrapolated value of J suffers a substantial 
uncertainty due to large finite-size corrections (that become even more sizeable for yet 
larger 7 values). 
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J 
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Figure 2. Asymptotic value of J as a function of the collision rate 7 for free BC 
and the same parameter values as in figure [TJ Data are plotted in log-log scales. The 
error bar for 7 = 5 has been estimated from a rough comparison among different 
extrapolation schemes. The dashed line is a power law fit of the first four data: its 
slope is -0.51. The solid line corresponds to the analytic solution for fixed BC (equation 
(20) of i). 

So far, we have not found any relevant difference between fixed and free BC. The 
heat flux scales in the same way in both cases and exhibits the same dependence on the 
collision rate. If the effect of the BC were restricted to a layer of finite width around the 
boundary, in the thermodynamic limit, the thermal resistance of a given chain would be 
independent of the type of thermal contact. In other words, we should expect J to be 
independent of the BC. However, this is not the case, as it can be inferred from figure 
[2J where we have also plotted the analytic curve for the fixed BC case (equation (20) in 
[6]). For free BC, the heat flux is approximately twice as that obtained for fixed BC. It 
is worth mentioning that the same effect was found in the simulations of FPU-/? chains 
[9], although with a slightly different value of the ratio (around 1.7 in that case). Since 
the flux is constant along the chain, this means that even deeply in the bulk, the system 
perceives the effect of the boundaries. In particular, from the knowledge of the local 
temperature profile and from the heat flux, one can in principle infer the type of BC. 
These results suggest that this is another way anomalous conduction manifests itself. 

The whole scenario is even more subtle than suggested by figure [21 In fact, for free 
BC, the leading term of the heat flux depends not only on 7 but also on the coupling 
strength A with the heat bath, while this is not so for fixed BC. We illustrate this in 
figure [31 where we plot the ratio 

j r =j{\ = l,N)/j(\ = l/4,N) (8) 
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Figure 3. Ratio of the energy fluxes as defined in equation JSJ, for free (diamonds) 
and fixed (circles) BC as a function of 1/ y/~N, and for 7 = 1 . The error on j r for free 
BC and N = 1600 is due to a not yet relaxed dynamics. 



where j(X, N) is the heat flux in a chain of length N and for a given value of A. It is 
not surprising to see that the coupling with the heat baths modifies the flux in chains 
of finite length. However, we see that for fixed BC, the effect of the coupling vanishes 
as j r converges to 1 (see the lower curve in figure [3]). On the contrary, for free BC, 
j r remains significantly different from 1. This suggests that fixed BC may lead to a 
kind of universal behaviour, namely the heat flux and also the temperature profile are 
independent of the details of the coupling with the heat baths. This is not the case 
for free BC. It would be interesting to check whether the same holds true in generic 
nonlinear chains. 

3.2. The temperature profile 

Another observable of interest is the temperature profile Tj = (pf). In figure H] we 
show the temperature profile for free BC and three different function of the 

"normalized" position along the chain y = 2i/N — 1, that varies in the interval [—1,1]. 
The data collapse is coherent with the scaling assumed in the continuum approach 
[6]. The shape of the profile is qualitatively similar to that obtained for fixed BC but, 
although here there are no (square-root) singularities at the boundaries of the chain (see 
equation (19) and (79) in [6]). Furthermore, we notice that the profile itself depends 
on both 7 and A. Evidence of such a dependence can be appreciated in the the inset 
of figure HI where the difference ST(y) between the profiles corresponding to A = 1 and 
1/4 is plotted in for three different system sizes. In fact, we see that ST(y) does not 
vanish in the thermodynamic limit. Moreover, the regions around the boundaries are 
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Figure 4. Temperature profile T(y) for 7 = 1 and A = 1 and free BC, for N = 100, 
200, 400 and 800. In the inset: the relative difference between the temperature profiles 
corresponding to A = 1 and A = 1/4, with 7=1 and for N = 400, 800, and 1600 
(dashed, solid, and dotted-dashed lines, respectively). 



affected by strong finite-size effects. In fact, one expects that 8T(—1) = ST(1) = 0, as 
the temperature necessarily converges, as N — > 00, to that of the attached heat bath. 

3.3. Other correlators 

In this section we analyse the behaviour of the different correlators ([3]), along the 
diagonal (i = j) and for generic values of x = (i — j)/VN^ We first analyse the 
case of fixed BC. 

At equilibrium, off-diagonal elements of the correlators ([3]) are zero. In the 
nonequilibrium steady state we have recently shown that the off-diagonal correlators 
v are of 0(1/ \/N) j6]. This is confirmed in figure [5^,, where we plot the lower diagonal 
of v^j+i, corresponding to v measured at a distance x = l/y/N from the diagonal, that 
is denoted by x = + . 

The potential energy profile y closely reproduces the kinetic energy profile (see 
also [5]). In order to appreciate its contribution, it is necessary to look at higher-order 
corrections. This can be done by introducing 

= - ^Vij ( 9 ) 

which measures the mismatch between kinetic and potential energy. Figure [5b shows 
that along the diagonal (see the lower set of curves), ip scales as 1/N everywhere except 
perhaps at the boundaries. From a physical point of view, this implies that everywhere 

|| The continuous coordinate x measures the distance of a given correlator from the diagonal. We have 
used the same notation in [6]. 
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Figure 5. Behaviour of the stationary covariances as a function of the coordinate y for 
fixed BC and different sizes: (a) First lower diagonal of v (b) diagonal (lower curves) 
and first subdiagonal (upper curves) of ip (c) diagonal of the symmetric component z + 
and (d) First subdiagonal of the antisymmetric component z . In all panels, dotted, 
dashed and solid lines refer to TV = 200, 400, and 800 respectively. The physical 
parameters are 7 = 1, A = 1. 



in the bulk, the system is locally at equilibrium (with 1/N finite-size deviations from the 
virial equality). The wild behaviour observed near the boundaries suggests the existence 
of nontrivial boundary layers. We will discuss this in detail in the next subsection. 
Analogously to v, ip exhibits a "discontinuity" when moving away from the diagonal. 
Indeed, in the upper set of curves of figure [5b we show that ip(x = + ,y) is of order 
1/N 3 / 2 . It is worthwhile remarking that this scaling holds only for A = 1. For A 7^ 1 
we have found that the off-diagonal terms of x/i are of order 1/N too. By recalling that 
we have here selected ou — 1 , it is reasonable to conjecture that the faster convergence 
of i/j observed for A = 1 is a manifestation of the thermal impedance matching on the 
boundary, theoretically predicted for A = uj (see [7j). 

Moreover, as shown in [7], it is convenient to distinguish between symmetric and 
anti-symmetric components of the correlators z with respect to x, 

z± = . (10) 

In figure \5jp, we plot the symmetric component which corresponds to the leading term 
of the heat flux. In fact, it scales as 1/y/N. The deviations from a perfectly flat 
shape reveal again the presence of nontrivial boundary layers. Along the diagonal the 
antisymmetric component z~ is zero by construction, while along the first subdiagonal, 
z~ scales as 1/N (see figure [5H). 

In figure [6] we show the behaviour of the correlators as a function of their distance 
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Figure 6. Off-diagonal behaviour of the stationary covariances for fixed BC and the 
same parameter values and same notations as in the previous figure: (a) v for y = —0.5; 
(b) i/j for y — 0; (c) z + for y = (d) z™ for y = 0. The physical parameters are 7 = 1 
and A = 1. For the sake of clarity, in panels a and b the diagonal values (x = 0) are 
omitted as they are of lower order. 

x from the diagonal. We have found that z~ and the derivative of z + along x are both 
discontinuous across the diagonal, in agreement with the theoretical analysis in [7j. All 
the results are independent of A except for the variable ip which, for A ^ 1, is constant 
away from the diagonal and of order 1/N . We would like to remark that this anomaly 
does not affect the theoretical analysis carried in [6] and |7J, as (for fixed BC) ip does 
not contribute to the leading behaviour of the temperature profile and of the heat flux. 

The very good overlap among the curves obtained for different system sizes confirms 
the scaling behaviour of the off-diagonal already seen in figure [51 Most important, the 
observed scaling corroborates the validity of the ansatz used in [6] and [7j. Summarizing, 
for fixed BC and far from the boundary we find that: 

- Along the diagonal, ip is 0(1/N). Off-diagonal, tp is 0(l/iV 3//2 ) for A = uj and 
0(1/ N) otherwise. 

- The correlator v is 0(1) along the diagonal and 0(1 /y/~N) off-diagonal. 

- The symmetric correlator z + is 0(l/\/~N) everywhere. 

- The antisymmetric correlator z~ is 0(1/ N). 

As a final remark, note that z + is the only variable that is continuous in x. This 
implies that the difference z + (0 + ,y) — z + (0,y), that we have denoted by 5z + in [7J, 
must necessary be an order e higher than its addenda, since its leading contribution is 
a derivative with respect to x. 
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Figure 7. Elements of the stationary covariances for free BC and different sizes: (a) 
z + (x = 0, y) (lower curves) and at z + (x = 0.7, y) (upper curves); (b) z + (x,y = 0); 
(c) ijj(x = 0, y) (d) z~(x,y = 0). In all panels, dotted, dashed and solid lines refer to 
N = 200, 400, and 800 respectively. 



We now turn our attention to the free BC. As it can be seen in figure [7J, the 
correlators scale with N in the same manner, irrespectively of the boundary conditions. 
We only notice the following qualitative differences: First, with free BC, the convergence 
at the boundaries is more effective than for fixed BC (compare figure [5b with figure [7b). 
Second, as a function of x, some additional oscillations of z can be seen only for fixed 
BC (figure [5b and figure [7b). Third, z + is larger for free BC, in agreement with the fact 
that in this case, the heat flux is about two times larger than for fixed BC. 

3.4- Behaviour at the chain edges 

The numerical discussion carried out in the previous subsection has revealed the 
existence of "boundary layers" in the vicinity of the contact points with the heat baths 
(y ~ ±1), where strong deviations from the expected scaling behaviour are clearly 
visible. Since in [7J, we have not attempted a theoretical analysis of the boundary 
layers, it is at least necessary to clarify their relevance, with reference to the numerical 
but otherwise exact solutions for the correlators. 

The variable that is mostly affected by the presence of such boundary layers is ip 
which even changes its scaling behaviour with N. This is shown in figure [8] for the case 
of fixed BC. In order to emphasize the scaling behaviour at the boundary, we subtract 
from if) the 1/N term (denoted by ^&) in the bulk that we know is constant (see [7J). 
For fixed BC and A = 1, ^ = 0, since the leading term is of order l/iV 3//2 , while for 
A = 0.25, ipb = 1/iV (with a few percent of uncertainty on the numerical constant.) The 
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Figure 8. Scaling behaviour of ipij close to the leftmost edge of the chain (y ~ — 1) 
for fixed BC and 7 = 1. Dotted, dashed and solid lines refer to N = 200, 400, and 800 
and A = 1. The three lower (upper) curves correspond to x = 0.7 (x = 0). Triangles 
correspond to N = 1600, A = 0.25, and x = 0. 



data collapse reveals that ip passes from values of order \/\N to values of higher order 
over a number of sites of order \/N. 

The existence of a boundary layer manifests itself in the values that different 
correlators assume at the boundaries (in the vicinity of y = ±1). At the level of the 
partial differential equations derived in [7], the BC (either free or fixed), lead to certain 
mathematical constraints among the correlators that must be satisfied for y — ±1. For 
instance, for free BC, we have found analytically that (see Eqs. (64), (65) and (66) of 

®) 

u 2 z + (x, -1) - Av(x, -1) = , and (u 2 - A 2 )v(x, -1) - u 2 tp(x, -1) = . (11) 

On the contrary, for fixed BC we have found that all correlators turn out to be zero at 
the boundaries. In figure [9] we have plotted the combined variables appearing in the 
l.h.s. of (fTTl) for uj and 7 unity. In panel a, we plot z + — Av as a function of y for A = 1 
(dotted and dotted-dashed curves) and A = 1/4 (solid and dashed curves). In both 
cases, the combined variable reaches zero at y — 1, thus confirming the findings in [7j. 
More importantly, for A = 1/4, one can infer that z + — Av will exhibit a discontinuity at 
y — — 1 in the limit — > 00. This is a direct consequence of the boundary layer which 
can be further seen in figure©, where we plot ip — (1 — A 2 )v for A = 1/4. Furthermore, 
we see that the second theoretical constraint ( TTTT) is also satisfied. This means that ip 
becomes of the same order as v. From a physical point of view this implies that in the 
boundary layer, i.e. at a short distance from the boundaries (of the order of y/N), local 
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Figure 9. The constraints (fTTj) among the correlators ip, z + and v are plotted as 
a function of y and w = 7 = 1. Dashed and solid curves correspond to N = 800 
and 1600, respectively, and A = 1/4. Dotted and dotted-dashed curves correspond to 
N = 800 and 1600, respectively, and A = 1. 



4. Relaxation to the stationary state 

Another interesting issue of nonequilibrium phenomena concerns the convergence 
towards the stationary state, with a particular reference to the time scales. In the 
context of our stochastic model, this question can be addressed by investigating time- 
dependent solutions of equations (j3J). As the evolution is linear, this can be done by 
determining the whole spectrum of the corresponding linear operator. 

For computational purposes it is actually convenient to recast the problem (J4]) in a 
more compact way, by a suitable "unfolding" of the elements of the matrices y, v and 
z in a linear array X. To minimize memory requirements, we take into account the fact 
that y and v are symmetric by construction and we consider only their independent 
entries. On the other hand, z is antisymmetric only in the stationary state . Therefore, 
at all finite times, all its elements must be considered. Altogether, X is composed of 
M = 2N 2 + N independent elements and the equations of motion for the correlators 
can be formally written as 

X = LX + X (12) 

where L is an M x M matrix and the vector Xq contains the source terms proportional 
to A. The matrix L is real but not symmetric and therefore, has M complex conjugate 
eigenvalues A = A# + iAj. Global stability of the stationary state requires all real parts 
A R to be non positive. 

The simplest approach consists in computing the eigenvalue spectra {Aj} of the 
matrix L with standard linear algebra algorithms. Their location in the complex plane 
is illustrated in figure [10] for three different values of 7. First we recall that, to our 
knowledge, this is the first nontrivial model where all time scales, from the microscopic 



Dynamics of anomalous heat transport: numerical analysis 



11 




Figure 10. Spectra of the matrix L for N = 20 and 7 = 0.5, 1, and 2 (panel a, b, and 
c, respectively) . For the sake of clarity we have used different scales along the vertical 
and horizontal axes with the exception of panel b, where we wish to draw the attention 
to the nearly circular symmetry. 

to macroscopic ones can be obtained at once. As expected, the whole spectrum lies on 
the negative A# semi-plane, confirming that the stationary state is stable. Secondly we 
note that the shape of the spectrum changes qualitatively upon varying the collision rate 
7. By increasing 7, the real part of the spectrum is shifted towards negative values. This 
is also not surprising as 7 quantifies the strength of the internal stochastic process and 
thereby of the corresponding relaxation processes. More interesting is the observation 
that the A^'s are distributed over an entire range of scales from 0(1) to very small 
ones. In the perspective of constructing a suitable hydrodynamic description (that is 
basically the goal of [7]), it is only the latter ones that matter. Unfortunately, we have 
not found a way to establish a direct connection between slow modes (those characterized 
by a small |A#|) and hydrodynamic modes, as this would require determining not 
only the eigenvalues, but also the eigenvectors. This task is numerically unfeasible, 
as the dimension of the space increases quadratically with N, and it is not even easy to 
determine the spectrum, let alone the eigenvectors. In practice, we have been able to 
determine the entire spectrum only up to N ~ 80. 

As far as we are concerned with the slowest relaxation processes, we can employ an 
alternative method akin to that used for the computation of the maximum Lyapunov 
exponent of a dynamical system. Indeed, for asymptotically long times 

X{t) = X(0)e Alt . (13) 

In order to estimate Ai, we integrated numerically the differential equations ( fl"2l) 
(actually it suffices to consider the homogeneous system X = LX since X is irrelevant) 
starting from random initial conditions with unit Euclidean norm, ||rY(0)|| = 1. For the 
sake of accuracy, we divided the time T tot of the whole run into n consecutive time 
intervals, each of length r, so that T tat /r = n. At the end of each time interval, we store 
the corresponding growth rate and renormalize the vector X to a unit norm. Finally, 
we determine Ai as the average 

1 n 

Ai = -E ln 11^(^)11- (14) 

T 1=1 
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Figure 11. Dependence of the leading eigenvalues on the chain length N, for 7=1 and 
fixed BC. Data is plotted in log-log scales. In panel (a), the maximum (real) exponent 
Ai (circles) is plotted together with a power-law best fit TV" 1,91 (dashed line). In 
panel (b), the real (circles) and imaginary (square) parts of the second eigenvalue A 2 
are plotted. The dashed lines correspond to the power-law fits N~ lm and _/V~ - 95 
respectively. 



As a result, we have been able to investigate systems of size up to N = 400. The 
numerical results plotted in figure [TTk show that for the considered parameter values, 
Ai is real and goes to zero with some power of N. A best fit suggests that Ai « N~ 2 . 
The same approach allows determining the corresponding (slowest) "mode" of the linear 
operator L. In figure [12] we plot the result obtained for N = 400. It looks very similar 
to the Fourier modes that we expect on the basis of the theoretical analysis carried out 
in [7j and the order of magnitude of the corresponding eigenvalue is in agreement with 
that analysis. 

For N < 80 we have been able to determine the entire spectrum. It turns out that 
the second and third eigenvalues are complex conjugate. In figure [TTb we plot their real 
and imaginary parts. The real part scales as 1/N 1 ' 91 , a value that, despite the limited 
amount of data, suggests again an asymptotic 1/N 2 behaviour. Instead, the imaginary 
part scales nearly as 1/N. Since this model is characterized by the presence of sound 
waves, we expect the imaginary part of A 2 to be connected to the periodicity due to the 
propagation of such waves. The period T of the oscillations can be written as 

T = 7T~T~ ~ ~ ' OS) 
(A 2 )/ c 

where c is the sound velocity (equal to ui in our arbitrary units), so that 
N(A 2 )j 

c = ~*r- (16) 

If we substitute for (A2)/ the value found numerically, we obtain c ~ 1, thus confirming 
our expectations. On the other hand, the 1/A 2 dependence of the real parts poses 
problems of consistency with the presence of an anomalous heat transport, as the 1/N 2 
dependence is expected to hold for normal heat conduction. In order get a better 
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Figure 12. Eigenvector of the linear operator L corresponding to Ai , for N = 400 
and 7=1. Starting from the highest panel on the left and going in clockwise sense, 
we plot the matrices yi,j, z 



?ij and z] - respectively. 



understanding of this, we have investigated the convergence of a specific observable, 
namely the flux J. More precisely, we have studied the relative deviation of the flux at 
time t from its asymptotic value (see equation (20) of [6]). This is shown in figure [TBI 
where, starting from an equilibrium state at temperature (T + + T„)/2 (so that J(0) = 0), 
5 r J(t) = 1 — is plotted as a function of time. We have considered two different 
definitions of instantaneous flux: i) the energy flux along the first bond (i.e. directly in 
contact with the heat bath), which corresponds to the dashed curves in figure [13l ii) the 
average flux (along the whole chain), which corresponds to the solid curves. Altogether, 
5 r J(t) is a measure of the deviation from the stationary state at time t. In order to 
compare the curves corresponding to three different sizes (N = 200, 400, and 800), 
the time variable has been suitably scaled (by a factor 8.95 for N = 200 and 2.97 for 
N = 400). The first part of the curves nicely overlap along a straight line, which signals 
an exponential convergence with a rate r](N) which, taking into account the temporal 
rescaling factors, is well reproduced by the law 

K") = + W < 17 > 

A best fit of the numerical results yields r/o ~ 4.4 in good agreement with the second 
eigenvalue of the operator theoretically derived in [7j, that is equal to 4.28 (see the 
spectrum plotted in figure 1 in [7], which has been obtained by setting all parameters 
equal to 1). The reason why the first eigenvalue does not play any role in our numerical 
study is that the corresponding eigenmode is not excited for our choice of the initial 
condition that is characterized by exactly the same average temperature as that of the 
asymptotic stationary state. 

The curves reported in figure [13] show that for any finite N there exists a crossover 
time beyond which a yet slower convergence sets in. By fitting the final slope, one can 




Figure 13. Relaxation of the energy flux to its stationary value for fixed BC, 7 = 1 
and A = 1. The relative amplitude of the heat flux (see the text for its definition) is 
plotted versus time for three different sizes (N = 200, 400, 800) and for two definitions 
of the flux. The dashed curves refer to the flux along the first bond of the chain; 
the solid curves refer to the spatially averaged flux. The time axis has been suitably 
rescaled to emphasize the scaling behaviour of the initial exponential regime. 



verify that the time scale of this last part of the convergence process is on the order of 
N 2 , in agreement with the previous spectral analysis. However, it is important to notice 
that this time increases with N and thereby corresponds to increasingly small scales 
(look at the vertical axis in figure [13]). Altogether, this means that the components of 
the initial state that lie along the slowest components become increasingly small upon 
increasing the system size N, until they vanish in the thermodynamic limit. One way 
to understand the unphysical character of these "super slow" modes is as follows. Any 
meaningful invariant measure is characterized by a set of correlators, but the converse is 
not true. Only the covariance matrices c that are positive semi-definite can correspond 
to physically meaningful state. For instance, we have verified that a sufficiently large 
perturbation along the eigenmode depicted in figure [12] leads to unphysical matrices. 



5. Discussion and conclusions 



The study of heat transport in a chain of N particles with nearest-neighbour coupling 
and conservative noise allows one to investigate both analytically [6] and numerically 
many subtle aspects of anomalous transport in one-dimensional systems. Since an 
analytical solution is available only for fixed BC, the free BC case can be investigated 
only by means of numerical methods. The comparison between the two cases shows 
that the physics of heat transport strongly depends on the choice of BC. For instance, 
as already observed in the FPU-/3 model [9] , the ratio between the heat fluxes measured 
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with free and fixed BC does not converge to 1 in the thermodynamic limit (N — > oo). 
Moreover, we find that for fixed BC the heat flux and the temperature profile are 
independent of the coupling strength with the thermal bath, while this does not hold 
for free BC. Nonetheless, the anomalous scaling of heat conductivity with the system size 
(k ~ iV 1 / 2 ) is found to be independent of the choice of BC. We have also investigated the 
convergence to the stationary state, both by determining the eigenvalues of the evolution 
operator of the covariance matrix, and by following the evolution of the average heat flux 
when starting away from the stationary state. The analysis reveals that over long time 
scales, the convergence is controlled by a rate r](N) which scales as N~ 3 ^ 2 . This means 
that if one wishes to extract reliable numerical data by performing direct numerical 
simulations, e.g., in a lattice of size N = 50000, it is necessary to evolve the system 
well above 10 time units. It should be kept in mind that similar limitations hold 
for deterministic nonlinear systems, even though they are, in general, characterized by 
slightly different exponents [9] . Finally, our analysis of the time-dependent solution has 
revealed a crossover from a typical fractional diffusion regime to a superslow relaxation. 
The crossover time is found to increase with the system size, suggesting that the latter 
regime becomes irrelevant in the thermodynamic limit. 
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